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1 Introduction and Summary of Results 

Reduction methods decrease the size and complexity of systems of kinetic 
equations. They are effective when a small number of variables can be sin- 
gled out as evolving on a "slow manifold" and the remaining (fast) variables 
somehow follow from the slow variables. In such cases, the system of kinetic 
equations can be reduced to a much smaller system for the evolution of only 
the slow variables, and the fast variables can be determined simply by table 
look-ups or by direct computation. Over the years, a large number of reduction 
methods have been proposed and implemented in computer codes; references 
can be found in our earlier article ^l] , and additional references are pH EH EOj • 

The focus of Ref. |Tl] was on the Intrinsic Low-Dimensional Manifold 
(ILDM) method due to Maas and Pope [d and an iterative method proposed 
by Fraser [H] and further developed by Roussel and Fraser [2E|. In this arti- 
cle, the focus is on the Computational Singular Perturbation (CSP) method 
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developed by Lam and Goussis [Zl El H21 H31 HH1 UHl CHI UHl 



A chemical kinetic equation is an ordinary differential equation (ODE), 



for a vector x of species concentrations; g is a smooth vector field, and t is 
time. Reduction methods are effective when the variables fall into two classes, 
fast and slow, as is the case when the Jacobian of the vector field has a spectral 
gap. For the analysis, it is convenient to identify the spectral gap with the 
inverse of a small parameter e, but we emphasize that this restriction is not 
necessary for the applicability of the CSP method. The characteristic time 
scales for the fast and slow species are given by the "fast" time t and the 
"slow" time r = et, respectively. We assume that the entries of x are ordered 
in such a way that the first m components evolve on the slow time scale and 
the remaining n components on the fast time scale. Then the vector field g 
has the form 



where I m and /„ are the identity matrices in R m and R ra , respectively, and 
the system (II. 1|) is a fast-slow system of ODEs. Both g 1 and g 2 may depend 
on e, but the entries of these vectors as well as their partial derivatives are all 
0(1) as e I 0, uniformly in x. 

Geometric singular perturbation theory (GSPT) jUHHI provides a natural 
framework for the analysis of fast-slow systems of ODEs. If such a system has 
a slow manifold, Aio, in the limit as e j and this manifold is normally 
hyperbolic, then GSPT identifies a (usually nonunique) slow manifold M. £ for 
e sufficiently small. GSPT also gives a complete geometric and analytical 
description of all solutions near M. £ , including how trajectories approach Ai £ . 
The goal of any reduction method is to find A4 £ , if it exists. 

Typically, the vector field g is written in a form suggested by chemi- 
cal kinetics, namely, as a weighted sum of the stoichiometric vectors, the 
weights being the associated reaction rates. But this representation is in no 
way unique. In fact, Eq. (jl.2|) shows an equivalent representation of g as a 
weighted sum of the standard basis vectors of R m+n , the weights being the 
coordinates egi, . . . , eg m , g m +i, . . . ,g m +n- The objective of the CSP method 
is to express g in yet another basis, one that is tuned to the dynamics of 
the system, where the fast and slow coordinates (amplitudes) evolve indepen- 
dently of each other. The CSP method achieves this objective constructively 
by successive approximation. Starting with a more or less arbitrary initial 
basis, one derives the evolution equations for the fast and slow amplitudes and 



dx 



(1.1) 




(1.2) 
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updates the basis iteratively in such a way that the evolution equations for the 
updated fast and slow amplitudes decouple to increasingly higher order in the 
small parameter e. Each iteration consists of two steps. The first step deals 
with the dependence of the fast amplitudes on the slow amplitudes, the second 
step with the dependence of the slow amplitudes on the fast amplitudes. 

After each iteration, one identifies the CSP manifold (CSPM) as the lo- 
cus of points where the then-current fast amplitudes vanish. The CSPM is 
an approximation to the slow manifold M. £ . The question is: How good is 
the approximation? In this paper, we analyze the general class of fast-slow 
systems of ODEs (|l.lj) - (|1.2j) and show (Theorem 13. lj) that the CSP method 
generates term by term the asymptotic expansion of the slow manifold Ai £ . 
After q iterations (q = 0,1,2,...), the asymptotic expansions of the CSPM 
and M. £ agree up to and including terms of 0(e q ); they differ in general at 
0(e q+1 ). Also, the qth application of the CSP algorithm leaves the terms at 
(9(1) through 0(s q ~ 1 ) invariant. (This observation is important because the 
lower-order terms have already been determined correctly in the preceding 
applications.) We illustrate Theorem 13 . 1 1 with an example from the Michaelis- 
Menten- Henri mechanism of enzyme kinetics [3J 122 E3 E2| • Similar results 
(for q = 1,2) have been obtained by Valorani, Goussis, and Najm [21] for a 
model equation due to Davis and Skodje [2j. 

Our proof proceeds via an intermediate result for a one-step CSP method. 
The one-step CSP method is the same as the full two-step CSP method but 
involves only the first step. It yields a sequence of slow manifolds, just like the 
full CSP method, whose asymptotic behavior as e j can be compared with 
that of the slow manifold Ai £ . The result ( Theorem 14. 1J) is that q applications 
of the one-step CSP algorithm yield an approximate slow manifold that agrees 
asymptotically with M. £ up to and including terms of 0(e q ). In other words, 
the one-step CSP method is as accurate as the full CSP method; and to prove 
the main result for the full CSP method, one needs only to show that the second 
step does not affect the lower-order terms in the asymptotic expansion of the 
CSPM. Although the second step of the CSP method does not play a role in 
the approximation of Ai £ , it does play a constructive role in the approximation 
of the dynamics near M. £l as we shall demonstrate in the special case of the 
Michaelis-Menten-Henri equations. 

In |TI], we showed that the ILDM method yields an approximate slow 
manifold that is asymptotically accurate up to and including terms of 0(e), 
with an error of 0(e 2 ) proportional to the curvature of Aio- The CSP method, 
on the other hand, can generate an approximate slow manifold that is asymp- 
totically accurate up to any order. The difference can be traced to two facts, 
namely, the choice of the fundamental operator governing the dynamics of 
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the system and the retention of the variation of the Jacobian over the man- 
ifold M.q. While the ILDM method is designed to transform the Jacobian 
of the vector field into triangular form (and often also into diagonal form), 
the CSP method is an algorithm to diagonalize the (nonlinear) Lie bracket 
involving the vector field to successively higher orders in e. The Jacobian is 
a linear approximation, so the ILDM method never gets beyond a linear ap- 
proximation. The variation of the Jacobian over M.$ introduces an extra term 
in the Lie bracket. By retaining it, the CSP method preserves the nonlinear 
character of the operator governing the dynamics of the system. 

This article is organized as follows. In Section we recall the Fenichel 
theory of GSPT and give the asymptotic expansion of the slow manifold A4 £ . 
In Section El we describe the full CSP method for fast-slow systems and state 
Theorem 13.11 The one-step CSP method is introduced in Section 0] The 
approximation result for the slow manifold is given in Theorem 14. 1| its proof 
occupies most of Section0]and uses two lemmas that are given in the Appendix. 
In Section we return to the full CSP method and prove Theorem 13.11 In 
Section |HJ we illustrate the CSP method and the results of this paper on 
a planar system of equations for the Michaelis-Menten-Henri mechanism of 
enzyme kinetics. Section [7| is devoted to a discussion of the relation between 
the CSP and ILDM methods. 



2 Fast-Slow Systems of ODEs 

Collecting the slow variables in a single (column) vector y and the fast variables 
in a (column) vector z, we rewrite Eqs. (jl.lj) - (jl.2j) as a fast-slow system, 

y' = egi(y,z,e), (2.1) 
z' = g 2 (y,z,e). (2.2) 

(A prime ' denotes differentiation with respect to t.) The long-term dynamics 
of this system are more naturally studied on the time scale of the slow variable 
r = et, where the system of Eqs. (I2.1j) - (j2.2|) assumes the form 

V = gi(y,z,e), (2.3) 
" = 92{y,z,e). (2.4) 

(A dot ' denotes differentiation with respect to r.) 

In the limit e { 0, Eq. ()2.4j) reduces formally to the algebraic equation 
#2(2/, z, 0) = 0. We assume that there exists a compact domain K 6 R m and 
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a smooth single-valued function ho on K such that 

g 2 (y,h (y),0) = 0, y e K. (2.5) 

Then the long-time dynamics of the system (|2.1|) - (|2.2j) are confined to the 
reduced slow manifold A4 , 

M = {(y, z) e R m+n : z = h Q (y), y e K}. (2.6) 

We assume, furthermore, that the real parts of the eigenvalues of the matrix 
D z g 2 (y,ho(y),0) are all negative, so M.q is asymptotically stable. Then the 
Fenichel theory which applies more generally to normally hyperbolic in- 
variant manifolds, guarantees that Aio persists as a slow manifold, so for all 
sufficiently small e there exists a slow manifold, M. £ , that is invariant under 
the dynamics of the system of Eqs. (j2.1|) - ()2.2J) . Moreover, M. £ has the same 
dimension as Aio an d lies near Aio, all nearby solutions relax exponentially 
fast to M. £ , and the long-term dynamics of the system ()2.1j) - (j2.2j) are governed 
by an equation on M. £ . The manifold M. £ is not unique; typically, there is 
a family of slow manifolds, all exponentially close (0(e~ c ^ £ ) for some c > 0). 
The following theorem is essentially a restatement of [10!, Theorem 2]. 

Theorem 2.1 For all sufficiently small e, there is a function h e such that the 
graph 

M £ = {(y,z):z = h £ (y),y£K} (2.7) 

is locally invariant under the dynamics of Eqs. \2.1)) - $Tty) . The function h £ 
admits an asymptotic expansion as e [ 0, 

he{y) = h (y) + ehiy) + e 2 h 2 (y) + ■■■, (2.8) 

and h £ G C r (K) for any finite r. The long-term dynamics of the system of 
Eqs. \2.1\) - En$l are governed by the equation 

V = 9i(vMv),e) (2-9) 
on M. £ , where ' = d/dr with r = et. 

The coefficients hi, h 2 , . . . are found from the invariance equation, 

g 2 {y, h £ {y),e) - e{Dh E ){y)gi{y, h £ {y),e) = 0, y E K, (2.10) 

in the following manner. (The invariance equation follows immediately from 
the chain rule, z' = Dh £ (y)y', and Eqs. ()2.1|) - ([2.2|h ) Each of the functions 
gi(- , h £ , e) and g 2 {- , h £ , e) admits a Taylor expansion near e = 0, 

CO oo 

gi{- ,h £ ,e) = ^2g hq e q , g 2 {- , h £ , e) = g 2>q e q , (2.11) 

q=0 q=0 
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with coefficients 

9i* = EE^r?(^^i)o Y,fa 1 ,.--,h ti ) + j{DZgi) 0i (2.12) 

k=0 j=l '3' \i\=q-k 

92, q = T,fa,.--,h ij ) + j{D*g 2 ) . (2.13) 

k=0 j=l \i\=q-k 

The notation ( • ) indicates that the quantity inside the parentheses is eval- 
uated on A4 — that is, at (y, h (y), 0). Note that (D J z D^g) is a multilinear 
operator, which maps a j-form to a vector. The inner sum in Eqs. (|2.12j) and 
(I2.13j) is taken over all multiindices i = (i x , . . . of j positive integers i\ 
through ij subject to the constraint |z| = i\ + ■ • • + ij = q — k. The expres- 
sions (|2.12|) and ()2.13|) hold for all q if it is understood that a sum is empty 
whenever its lower bound exceeds its upper bound. Substituting the expan- 
sions (J2.12|) and (J2.13j) into the invariance equation, Eq. (|2.1(J|) . and setting 
the coefficient of e q equal to zero, we obtain an infinite set of equations, 

<?-i 

92, g - ^(Dh^-i-t = 0, q = 0, 1, . . . . (2.14) 

1=0 

The first few equations are 

92 fi = 0, (2.15) 

{D z g 2 ) hi + (D e g 2 ) - {Dh )g 1>0 = 0, (2.16) 
(D z g 2 ) h 2 + \{D 2 z g 2 ) Q (h u h) + {D s D s g 2 ) Q hx + \{D 2 £ g 2 ) 

- (Dh 1 )g lfi - (Dh ) {{D z g x ) h x - {D £ g x ) Q ) = 0. (2.17) 

Equation (j2.15J) is satisfied identically, Eq. ()2.16|) yields the coefficient h x , 
Eq. ()2.17|) the coefficient h 2 , and so on. 

Remark. The assumption that the chemical species can be divided into fast 
and slow species, as in Eqs. (I2.1|) - (j2.2|) . is made for convenience. Our analysis 
can also be applied to general chemical systems where each species may be 
involved in both fast and slow reactions and for which there is a slow manifold. 



3 The CSP Method for Fast-Slow Systems 

In Eq. (I1.2|) . the vector field g is represented in terms of the standard orthonor- 
mal basis. It is useful to examine the representation of g in terms of other 
bases, especially bases whose entries depend on x. 
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Let A be an (m + n) X (m + n) matrix whose entries may depend on x 
and whose columns form a basis for the space R m+ri for each x. The vector 
field g may be expressed in terms of this (variable) basis A as 

9 = Af, (3.1) 

where / is the vector of the coordinates (amplitudes) of g. (When the columns 
of A are the stoichiometric vectors, the amplitudes are the reaction rates.) 
Since A is invert ible, B = A -1 , and 

/ = Bg. (3.2) 
The amplitudes can be split into two classes, / = I f2 , where f 1 is an 



. f . 

n- vector representing the fast amplitudes and f 2 an m-vector representing the 
slow amplitudes. The splitting suggests that we consider a decomposition of 
A, namely, A = (Ai,A 2 ), where A\ is (m + n) x n and A 2 is (m + n) x in, 



and a corresponding decomposition of B, namely, B = ^ 2 j , where B 

is n x (m + n), and B 2 is m x (m + n). Thus, f 1 = B 1 g and f 2 = B 2 g. 
Also, the identity AB = I on R m+n implies that A X B X + A 2 5 2 = J, while the 
identity BA = I on R m+n implies that B 1 A 1 = I n and 5 2 Ai = on R n and 
B 2 A 2 = I m and B l A 2 = on R m . 



The fast and slow amplitudes evolve in time. Differentiating Eq. ()3.2j) 
along solutions of the system (jl.ljl . we obtain 

— = B— + — ^ = B[Dg)g + — g 
c?t dt dt dt ' 

where Dg is the Jacobian of g. Hence, / satisfies the nonlinear ODE 

| = A/, (3.3) 

where A, the generator of the dynamics for the amplitudes, is given by 

d R 

A = B(Dg)A + —A. (3.4) 

Since BA = I and / is time invariant, A, B, and their time derivatives satisfy 
the identity 

(dB/dt)A + B(dA/dt) = (3.5) 
at all times. Hence, the definition (|3.4j) is equivalent to 

d A 

A = B(Dg)A — B—, (3.6) 
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where dA/dt = (DA)g. For completeness, we note that the identity (|3.5J) 
implies that {{DB)Af)A + B((DA)Af) = 0. 

In general, the operator A is not diagonal, and the equations governing 
the evolution of f 1 and f 2 are coupled. An ideal basis A is one in which A 
is block-diagonalized, so that the ODEs for f 1 and f 2 decouple. The CSP 
method approaches this ideal by successive refinements of the basis matrices 
A and B. The algorithm starts from a constant matrix A^°\ 

*«> = K> 4°>) = ( |! J I ■ (3.7) 

Here, A^ is an m x n matrix, A^ annxm matrix, and the off-diagonal blocks 
A^ and A$ are full-rank square matrices of order m and n, respectively. A 
common choice is A$ = 0, so every column vector of A± lies in the fast 
subspace. We follow this convention and assume, henceforth, that A^ = 0, 

A0» = (A?\ 4°») = ( ° $ ) • (3.* 

^21 ^22 

A more general choice of y4^ ^ is discussed below, after Theorem 13.11 The 
inverse of A<°> is 



5(0) = (^ (0) )- 



5 



1 

(0) 



■(A^)-^(Af^ 



*12 



R21 

5 (o) 


Rl2 




(4?: 










(3.9) 



The algorithm then proceeds iteratively. For q = 0, 1, ... , one first defines 
the matrix A( 9 ) in accordance with Eq. ([3. tip . 

A (9) = - = f S !l ) , (3-10) 

and matrices and L( g ), 



■(""Vo o "<«» - UtoVNi)" 1 °. 

Then one updates A^ and Bi q \ according to the formulas 

AW = A^(I-U {g) )(I + L {q) ), (3.12) 
B (q+1) = (I-L {q) )(I + U {q) )B {q) , (3.13) 
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and returns to Eq. f)3.10j) for the next iteration. 



At each iteration, one imposes the CSP condition, 

B\ q) g = U, g = 0,l,..., (3.14) 

to identify those points where the fast reaction rates vanish with respect to 
the then-current basis. For q = 0, is constant and given by Eq. for 
q — 1,2, ... , the CSP condition takes the form 

B \ q ){y^(q-i){y^)^)9{y,z,e) = ^, g = i,2, .... (3.15) 

If, for any q, the CSP condition is satisfied by a function z = ip(q)(y,£), then 
= {(y,z):z = ij {q) {y, e), y e K}, q = 0, 1, . . . (3.16) 
is defined as the CSP manifold (CSPM) of order q. 

Theorem 3.1 The CSP manifold K,^ agrees asymptotically with A4 £ up to 
and including terms of 0(e q ) for q = 0, 1, . . . , 

i 

iJ iq) (-,e) = J2^h, + 0^ +1 ), elO. (3.17) 

3=0 

Our proof of Theorem 13.11 proceeds via an intermediate result, which is 
of independent interest. We introduce a "truncated" CSP method, where we 
apply, at each iteration, only the first of the two steps of the full CSP method 
and skip the second step. This one-step CSP method reduces the matrix A 
to lower block-triangular form. We show that, after q iterations, the one-step 
CSP method generates a manifold ic!f \ whose asymptotic expansion agrees 
with that of Ai e up to and including terms of 0(e q ) (Theorem I4.1J1 . In other 
words, the one-step CSP method is as accurate as the full CSP method is 
claimed to be in Theorem 13.11 We then return to the full CSP method and 
carry out an asymptotic analysis of the modifications introduced by the second 
step. This second step reduces A further to block-diagonal form. We show 
that, at the qth iteration, the second step affects only terms of 0(e q+1 ) and 
higher. Hence, K!f^ approximates Ai £ as accurately as K.i 9 \ and Theorem 13. II 
follows. 

Theorem 13.11 extends readily to the case where the eigenvectors of the 
Jacobian Dg are used, instead of the stoichiometric vectors, to form the initial 
basis A(°\ In that case, the slow subspace of the leading-order Jacobian coin- 
cides with the tangent space T V M.§ at any point p G A4q, so the columns of 
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A 2 are tangent to M.q to leading order. In turn, this implies that the rows of 
.BL (p) span the orthogonal complement of the tangent space, also to leading 
order. As a result, the initial CSPM, the solution of Bh,g = 0, coincides with 
M. £ up to and including terms of 0(e), which is one order higher than is the 
case when A^ ' is given by Eq. (|3.8|) . Moreover, for each q = 1, 2, ... , the proof 
of Theorem 14.11 generalizes directly to this case. The asymptotic expansion of 
if)r q \ coincides with that of h e up to and including terms of 0(e q+1 ), which is 
one order higher than is the case when is given by Eq. (|3.8|) . 



Remark. Lam and Goussis, in their presentation ^2] of the CSP method, 
perform the update 1)3.12)1 and (J3.13)) in two steps. The first step corresponds 
to the postmultiplication of A^ with I—Ur g \ and premultiplication of B^ with 
I + U( q ), the second step to the subsequent postmultiplication of A™ (I — U^) 
with J + L( g ) and premultiplication of (I + U^)B^ with I — L^. The nonzero 
entries of U^ q ) and are chosen so that A is block-diagonalized to successively 
higher order in e. 



Remark. The definition ()3.6j) implies that A is the product of B with the 
Lie bracket of A (considered column by column) and g, 

A = B [A,g] = B([A. fl ,g], . . . , [A. tm+n ,g]). (3.18) 

The Lie bracket of two vector fields a and g is [a, g] = (Dg)a — (Da)g |21j . 



Remark. It is useful to state how A transforms to understand its properties 
as an operator. If A = AC and B = C~ X B, where C is an invertible square 
matrix representing a coordinate transformation in R m+n , then 

A = B(Dg)A-B^ = C- 1 B(Dg)AC-C- 1 B^^- 

dt dt 

= c- 1 b(d 9 )ac-c- 1 b(^c + a^-) 

\ dt dt J 

dC 

= C^AC -C- 1 ^, (3.19) 
at 

where dC/dt = (DC)g. The presence of the term C~ x dCjdt in Eq. ()3.19J) 
shows that A and A are not similar unless C is constant. 
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4 The One-Step CSP Method 



The goal of the one-step CSP method is to reduce the matrix A to lower 
block-triangular form — that is, to push the matrix A 12 to increasingly higher 
order in e. The method is identical to the full CSP method except for the 
updating of the matrices A and B. One starts from the same bases, = A^ 
and B(o\ = -B(o), and, instead of Eqs. (j3.12|) and (j3.13|) . uses the one-step 
expressions 

= A^(I-U {q) ), (4.1) 
B (q+l) = (I + U {q) )B {qh (4.2) 

where the matrix U( q ) is defined as in Eq. (|3.11|) with A replaced by A. (A 
tilde ~ distinguishes a quantity from its counterpart in the full CSP method.) 

The update rule for A follows immediately from Eq. ()3.19|) . 

Afe+i) = (J + Efa) A W (I - U {q) ) + (I + U( q) )^. (4.3) 

(Note that the identities v4 (0) = A (0) and £? (0 ) = B {0) imply that A (0 ) = A (0 ).) 
The matrix and its time derivative have the same block structure; only 
the upper right block is nonzero, so U^dUt q \/dt = 0, and Eq. ()4.3|) reduces to 

A (9+1) = (/ + U (q) )~A (q) (I - U (q) ) + (4.4) 
In terms of the constituent blocks, we have 

( 4 - 5 ) 



A 11 


-1) 




A 12 
A (<H 


-1) 


= ^(<?) A (J " 


121 
A (<H 


-1) 


_ 121 
- A ( g ); 


X 22 
A (</H 


-1) 


A 22 a! 
- A (g) - A ( 



[21 r7 i jjM 



21 JT 

(q) U (<l)i 



(4.6) 

(4.7) 
(4.8) 



where we have used Eq. (|3.11j) to simplify Eq. (|4.6|) . Note that we freely use 
U( q ) to denote both the full update matrix and its restriction to the subspace 
R m ; the latter is represented by the matrix (A/\) _1 AR. The appropriate 
interpretation is clear from the context. 

The one-step CSP method generates a sequence of manifolds, 

iCP = {{y,z):z = $ iq) {y,e),yeK}, q = 0, 1, . . . , (4.9) 
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just like the full CSP method; cf. Eq. ()3.16|) . The functions ip(q) are defined 
by the conditions 

B\ q) g = 0, g = 0,l,..., (4.10) 
where Bl q , is obtained from Eq. (|4.2|) . 

Theorem 4.1 The manifold K,^ agrees asymptotically with Ai £ up to and 
including terms of 0(e q ) for q — 0,1, ... , 

q 

h q) (-,e) = Y / e j h j + 0(e 9+1 ), e 10. (4.11) 
The proof of the theorem is by induction on q. 



4.1 The Induction Hypothesis 

The central idea of the proof of Theorem 14.11 is to express the CSP condi- 
tion (|4.1()jl in a form that resembles that of the invariance equation (|2.1()jl and 
then to derive the conditions under which the left and right members of the 
two equations are the same at each order. 

We begin by expressing the quantities A^ q+l \ Br g+ ^, and A( g+1 ) in terms 
of the original quantities A^°\ Br \, and A( ). Applying the definition (|4.1j) 
recursively, we find 

A^ = A^f[(I-U U) ). 

Since each Uu\ is nilpotent, it follows that 

A {q+l) = A^(I-P (q) ), (4.12) 



where 



Similarly, 



4)=E^)=( I E? - (A «'» r A «'» )■ (4.13) 

3=0 



B (q+1) = (I + P {q) )B {0) . (4.14) 

Substituting Eqs. (|4.12j) and (|4.14|) into the transformation formula ([3. 19)1 . 
and recalling that A( ) = A( ) and P^dP^/dt = 0, we find 

A (g+1) = (J + 4))A(o)(/ - P iq) ) + (4.15) 
12 



We use these expressions to rewrite Eq. (14. 10(1 . Since B?£ = 0, the equation 
becomes 



B$g 2 + e P^&fa + B 



21 



)11 

(0) 



or, since B}1 = 



'(0) 



92 + eA 21 P( 5 _ 1 )5( 1 ) + 5 (0) 



01 = 



<?i = 0. 



(4.16) 



The last equation has the same form as the invariance equation (|2.10j) . The 
solution of Eq. (|2.1U|) is z = h e (y), which defines A4 e , while the solution of 
Eq. (|4.1(j|) is z = ^}r q )(y,e), which defines JCi Q \ 

We analyze the CSP condition (|4.16|) order by order, up to and includ- 
ing the terms of 0(e q ). We recall that the components of the vector field 
g(y,z,e) are evaluated at z — ipr g \(y,e), the matrix P( q -i) is evaluated at 
z = i>( q ^i)(y, e), and the blocks of and Br \ are constant. Substituting the 
asymptotic expansion of ip(q), 



(4.17) 



into Eq. (|4.16|) and setting the coefficients of l,e, ... ,e 9 equal to zero, we 
obtain a set of equations, 



■P(«-1,0) 5 (0) + S (0) S'l.j-l + X] A 2lP(q-l,i) B (0)9l,j-£-l = °> ( 4 - 18 



£=1 



for j = 0, 1, . . . , q. Here, -P( g -i/) is the coefficient of the C(e £ ) term in the 
asymptotic expansion of P( q -i)- 

Equation ()4.18|) defines V'(gj) f° r ^" = 0, 1, . . . , g. The leading-order (j = 0) 
equation in the system (J4.18|) is the same for all q, 



92(y,4>( q ,o)(y),o) = o, <? = o, 1, 



(4.19) 



This is also the equation defining h . Its solution need not be unique, but we 
can identify each ij)(q,o) with h , 



ijj(q,o){y) =h (y), g = 0,l,.... 



(4.20) 



Then also i/j( q )(- > 0) — for g = 0, 1, . . . , so to leading order each manifold 
K$ coincides with A4 . 
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We wish to show that i>{ q ,j) = hj also for j — 1, 2, . . . , q. To this end, we 
compare Eqs. ()2.14|) and (j4.18|) . For a fixed j, the two equations match if 



A 



(0) 
21 



P {q - m B^ + B 11 



(0) 



AO) p R 21 

A 2i Mg-i^)-D(o) 



-Dh f , 



(4.21) 
(4.22) 



Conversely, if Eqs. (jQH) and hold, then = /i,,-. Notice that 

Eqs. f!4.21|) and ()4.22|) are independent of j; hence, they are nested, in the 
sense that, when j is increased by one, the equations for lower values of j 
remain the same. Thus, it suffices to prove Eqs. (|4.21j) and 1)4.22)1 for j = q. 
The proof is by induction on q, where the induction hypothesis is 



.4 



(0) 
21 



P {q -i ) {-,'ip {q -i ) ,e)B^ ) + B 



)ii 

'(0) 



9-1 

: -^eWhj + Oie 11 

j=0 

1 

4)(-,e) = "£eih j + 0(e« +1 ). 

j=0 



(4.23) 
(4.24) 

(4.25) 



The validity of these equations for q = 1 is shown in Section l4~2l The induction 
step is carried out in Section l4~3l 



4.2 Proof of Theorem 14.11 for q = 1 

We fix q = 1 and consider the 0(e) terms of Eq. ()4.16j) . 



(0) 



(£^2)0^(1,1) + (D e g 2 )o + A\ x 



-P(o,o)-B( 2 o) + B 1 ^ 



0i,o = 0. 



(4.26) 



The first and second terms in this equation are exactly the same as those in 
the equation for hi, see (I2.16)) . Therefore, we need only to show that the third 
term equals —(Dho)gi$ in order to prove the theorem for q — 1. 



According to the definitions ()4.13)) and ()3. 1 1|) with q = 0, we have 



Pi 



(0) 



(0) 



(A 



11 S-l A 12 

(0)J A (0) 



(A 



11 \-1a12 
(0)J A (0)> 



(4.27) 



where A( ) = B^(Dg)A^°\ according to the definition in Eq. ()3.1()j) . Now, 
A(o) admits an asymptotic expansion, A( ) = YljLo £ ''M d)i anc ^ eacn °f the 
coefficient matrices A(oj) consists of four blocks, 



A 11 

A (o,i) 



A 



(0) 

21 5 



(4.28) 
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A 12 
A (oj) 



A 21 

A (oj) 



B 



12 

(0) 



(D y g 2 ) 3 A { $ + (D z g 2 ) 3 A 



(0)' 
22 



+ B 



ii 

(o) 



(o) 

22 



^A^);-i4°i\ 



A 



22 

(OJ) 



5 



21 

(0) 



(0)' 
22 



(4.29) 
(4.30) 
(4.31) 



The notation ( • )j indicates the jth term in the asymptotic expansion of the 
quantity inside the parentheses, and it is understood that such a term is absent 
if the subscript is negative. 

A direct evaluation shows that the blocks Am^ and A} 2 m are nonzero. 
Therefore, A}L and A/£ are both 0(1), and 



(0,0) 



1 A 12 



A 



B 



12 



L (o,o) 



(D z92 )- 1 {D vga ) Q A<g +4? • (4-32) 



70,0) - ^(o,o) - (A( 0) o); iS -(o,o) - ^(o) 
Here, all the quantities are evaluated on Aio, where the identity 

(D z g 2 ) 1 (D y g 2 ) = -Dh 
holds. Hence, Eq. fl4.32|) implies 



|(0) 



P, 



(0,0) 



(0,0) 



(0) 



(4.33) 



(4.34) 



Finally, substituting this expression for P(o,o) hito Eq. ()4.26|) and using the 
identity A^B 11 



io) - -^aa-Bjo), we obtain 



(-0*02)0^(1,1) + (D s 92)o ~ (Dh )gi fi = 0. 



(4.35) 



This equation for ^(1,1) is the same as Eq. ()2.16|) for hi; hence, ^(1,1) = hi and 
V'(i) = ^0 + ehi + 0(e 2 ). This proves the theorem for q = 1. 

4.3 Proof of Theorem 14.11 for q = 2, 3, . . . 

We prove that Eqs. ()4.23|) - (j4.25|) hold for q + 1, assuming that they hold for 
0, 1, . . . ,q. By our discussion of Eqs. (|4~2"T]) and Eq. (jOS)) follows im- 

mediately from Eq. (|4.24|) . so we need only to consider Eqs. ()4.23|) and ()4.24|) . 



4.3.1 Establishing Eq. fflTZSh 

We first consider Eq. (|4.23jl . The induction hypothesis gives the estimate 
U {i) {- ,^ (i) ,e) = 0(e*) for i = 0,1, . . . ,q - 1. Also, ^ } = + G(e i+1 ) for 
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% = 0,1, . . . ,q-l. Hence, U^(- ,V>( g ),e) = U^(- ,ip(i),e) + 0(e l+1 ), from which 
it follows that 

U®(- , 4), e) = 0(e l ), i = 0, 1, . . . , q - 1. (4.36) 
In particular, U(p)(- ,ip( q ),e) = O(l), so 

1=0 

This asymptotic estimate can be used to derive asymptotic expansions of the 
blocks of A( g ). We begin with A}\. From Eq. ()4.15|) . we have 

AJJ = Ag ) + P (ff _ 1) Ag ) . (4.37) 
Since A£k = 0(e) by Eq. (|4.30j) . we see immediately that 

\ 9 ) = ^(0,0) 

Next, we examine the block A^. From Eq. (|4.6|) . we have 



AJJ = A" +0(e). (4.38) 



A (g) = f/ (9-l) A (g-l) - f/ (9-l) A (g-l) f/ (9-l) + 



First, [/(,_!)(• ,4), e) = O^ 1 ) by Eq. (JO^. Also, Ag_ 1} = Ag } = 0( e ) on 
by Eqs. (jHU) and (jOU]). Moreover, Ag_ 1} = Ag - AgjP^a) = C?(e) 
by Eqs. ()4.15|) and (|4.31|) . Finally, by applying Lemma lA.21 with V = 
we find that dU( q -i)/dt is 0(e q ). Putting these estimates together, we obtain 
the estimate 

A S = ^ A L) + (^ +1 )> ( 4 - 39 ) 

where we grouped all of the 0(e q ) terms into e q A}^ q y By combining the 
definition (|3.11|) with Eqs. ()4.38|) and ()4.39|) . we derive the desired estimate, 
U iq) = ( A ;ir lA S = 0(e q ). 



Remark. While the estimates of A?i and A?^ are not needed here, they will 
be needed in Section B First, A?J = A£k = O(e) on by Eqs. gJ5|) 

and g3DJ. Then, Af g) = Ag - A^P^ by Eq. (JHHJ). Now, Ag = 0( e ) 
by Eq. (|4.31j) . and thus the discussion for the size of A^ also yields that 
A f«i = ®( £ )- Putting the estimates of this section together, we obtain 



An(. *ne)=\ A (io) + ^) -% 2 9 ) + ^ +1 ) ) (440) 
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4.3.2 Establishing Eq. (JHHD 



Next, we consider Eq. ()4.24|) . The induction hypothesis gives the estimate 

4? [P® (" , hh £ ) B (o) + B (o)} = ~ E;=o £jDh J + for i = 0, 1, . . . , g - 1. 

Our goal is to show that this equation also holds for i = q. We first show that 
the terms up to and including 0{e q ~ 1 ) in both members of the equation agree 
for i = q. Then we analyze the terms of 0(e Q ). 



By the induction hypothesis, we have the asymptotic expansion 



.4 



(0) 
21 



P iq -i ) (-,i> iq -x),e)B^ ) + B 



; (o) 



9-1 
3=0 



(4.41) 



Also by the induction hypothesis, ip{q) — i>(q-i) + @(. £9 )- Hence, 



.4 



(o) 



21 



)11 

'(0) 



p { q -i)(->A*)> e ) B $) + B 

The definition (j4.13J) of P( 9 ) yields the update formula 



9-1 

3=0 



(4.42) 



(<?) 



P(q-l) +U (q ). 



(4.43) 



We already showed that U< q ){- , ip( q ), e) = 0(e q ), so Eq. ()4.43|) implies that 
the asymptotic expansions of P( q )(- , V'(g), £ ) an d P( g _i)(- , i[)( q ), e) agree up to 
and including terms of 0(e q ~ l ). The same, then, holds for the asymptotic 
expansions of A$ [P {q) (• , $ (ff) , e)B^ ) +B^ ) ] and A$ [P( q -i) (• , tf( q ), e)Bgj+.Bj|j]. 
Therefore, 



,(0) 
^21 



In other words, 



P (9) (.,V(,),e)S^ + S 



'(0) 



9-1 

i=o 



(4.44) 



.4 



(0) 
21 



-Pfeo)^ 1 ) + -B("o) 

4 (0) A -n R 21 
^21 - r (g,j)- D (o) 



—Dhj, for all J = 1, 



which establishes Eq. (|4.24j) for all terms up to and including 0(e 



(4.45) 
(4.46) 



It remains to show that the terms of 0(e q ) in both members of Eq. (|4.24jl 
agree, that is, 



A {0) P, ,R 21 



-Dh q . 



(4.47) 
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We achieve this by deriving an explicit formula for A^i P( q , q )B 2 ^ and com- 
paring it to that for Dh q , which is given in the Appendix (Lemma lA.lj) . 
We proceed in two steps. In step one, we express P(q, q )B 2 ^ in terms of 
P( q -i,o), ■ ■ ■ ,P(q-i, q -i)- Then, in step two, we obtain the explicit formula for 



(9-1,0) > 

A^i P( q , q )B^ in terms of the vector field and of Dhi, i — 0, 1, 



q-l. 



Step 1. Recall the update formula ()4.43|) . P( g ) = P( g -i) + U( q )- Using 
the definition (|3.11|) of and the explicit formula (|4.15j) for A( g ), we can 
express in terms of A( ) and P( q -i)- In particular, Eq. (|4.4U|) implies that 
U(, 



'(9,9) = ( A (o,o))~ A fe9)- Also > E( l- gives 



A (g) = A (0) - A (0) P (9-1) + P (9-1) A (0) ~ P (9-1) A (0)^(9-1) + 

It follows that 

U {q,q) = ( A (0,0)) _1 



dP, 



(9-1) 



(It 



A (0, 9 ) - ( A (0) P (9-1) 



(^(9-i)A?o 2 )) ( 



-(^(9-1)^0)^(9-1)) + 



dP, 



(9-1) 



dt 



(4.48) 



(4.49) 



where we recall the notational convention that ( ■ ) q stands for the coefficient 
of the 0(e q ) term in the asymptotic expansion of the quantity in parentheses. 
Using Lemma fA. 21 with V = P( q -i) and the fact that A 2 ^ and A^ ^ are both 
zero, we rewrite Eq. (|4.49J) as 



(9,9) 



(Afto))" 1 



Jl+{J% — A (0,0)-P(9-l,9)/ 



J% + J A + Jr 



where 



9-1 



9-1 



Jx = A 12 



(0,9)' 



~y] A ?0,q-O-P(g-l^)» 



J. 



-P(9-l,^)A(0, g - 



£=0 



q— 1 q—X—i 
i=0 j=0 



£=0 
9-1 



21 p 

(O.g-i-^M?-!,*)) 



£=0 



(9-l.Q 

dy 



9i,a-i- 



(4.50) 



(4.51) 



Substituting the expression ()4.50|) into the update formula (|4.43jl for P( 3 ), we 
find 



^21 P{q,q)B (0) 



^(Afjo))" 1 [ J i + J 2 + J 3 + ^4 + J 5 ] B 21 



(0)' 



(4.52) 



Step 2. We rewrite the terms Ji, . . . , J5 by means of the induction hypothesis 
and the explicit formulas ()4.28|) - ()4.31|) for the blocks of A( ). 
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Equation (|4.28jl and the identity A^B}^ = I n imply that 

(^lo))' 1 = ((DM' 1 Afl (4.53) 

Here, (D z g 2 )o stands for the leading order term in the asymptotic expansion of 
(D z g 2 )(- ,i/)(q),s). Since ip{q) an d h e agree up to and including 0(e q ) terms by 
assumption, the asymptotic expansions of (D z g 2 )(- , i>{ q ), s) and (D z g 2 )(- , h e , e) 
also agree up to and including 0(e q ) terms. For the remainder of this section, 
it does not matter whether quantities are evaluated on K.i or on A4 £ , since 
only the coefficients of e q or lower appear in our formulas. Accordingly, we 
make no distinction between the asymptotic expansions of a quantity evaluated 
on the two manifolds. 



Using Eq. g^Sj) and the identities B^ } = (A^) -1 , Bgj = (^f?) -1 , and 
flg = -B^A^B^, we find 



4? JiB$ = (D y g 2 ) q + {D z g 2 ) q A^B^ 

- A { SB^(D y9l ) q ^ - A^B^g^A^B^. (4.54) 

Next, substituting for A 2 } P^i^B^k from the induction hypothesis (|4.24|h we 
obtain 

3-1 

A^J 2 B^ = Y,{D z g 2 ) q -,Dh t - {D z g 2 ) q A^Bf ) 
e=o 

9-1 

-Y J A 2 Q 2 ) Bll ) {D z g 1 ) q ^Dhi + Af 2 ) B% ) (D z g 1 ) q .- 1 Af 2 ) B% ) . (4.55) 
e=o 

Then, using Eq. (j4.31|) and the assumptions of the lemma, we find 

3-1 a-i 
A^JsBtu) = -Y,Dht{D y g x ) q ^ -^^(D^^A^B^ 

+ A^B^iDyg,)^ + A^B^(D zgi ) q ^A^B^. (4.56) 
In the same vein, we use the induction hypothesis on J 4 , 

q— 1 q— 1— i q— 1 

A§}JaB^= Dh^D^^^Dh^Y.^Bt^D^q-i^DK 

i=0 j=0 i=0 

9-1 

+ Dhj(D z gx) q -i-jA 22 B(ty - A^Bf^D^^A^B^. (4.57) 

3=0 
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The terms in Eq. 1)4.52)1 containing A^ 2 sum to zero, which may be seen as 
follows. The second and fourth terms in ()4.54j) cancel against the second and 
fourth terms in ()4.55)) : the third term in ()4.54|) cancels against the third term in 
(14.56)1 ; the third term in (J4.55)) cancels against the second term in (j4.57j) ; and 
the second and fourth terms in (j4.56j) cancel against the third and fourth terms 
in ()4.57j) . These cancellations were to be expected because the approximation 
should be independent of the choice of A^°> . 



Carrying out the same type of calculation as above, we obtain 



a(°) T R21 
/i 21 ^5-0(0) 



9-1 

^2D 2 h e g liq ^ 



(4.58) 



where we have used the symmetry of the bilinear form D 2 hi. 

Equations ()4.53j) - ()4.58)) . together with the observed cancellations, yield 
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((^2)0 



-1 



j(0) p D 

A 21 r {q,q) n {0) 

q-1 q-1 

{D y g 2 ) q + S ^2{D z g 2 ) q -tDhz - D 2 h e gi^^ 

q—l q—l q — l—i 

z 

9i) q -i-i-jDhi 

1=0 i=0 j=0 



(4.59) 



A term-by-term comparison with the expression for —Dh q given in the Ap- 
pendix, Eq. (|A.3j) . shows that A^ P^^A 2 ^ = —Dh q . Thus, the proof of 
Theorem 14. II is complete. 



Remark. In general, the error term is nontrivial, as can already be seen at 
q = 0. The equation determining ^(0,1) * s 

{D z g 2 )ai> m + (D £ g 2 ) Q - A^B^ )9l>0 = 0. (4.60) 

This equation is not the same as Eq. (j2.16|) . which determines h\. Where 
Eq. (|2~TH|) has the term Dh , Eq. (ICTl) has the term A$B$ When the 

slow manifold is nonlinear, Dho depends on y, whereas A 22 B?^ is a constant 
matrix. Therefore, in general ^(0,1) 7^ hi, an d the strongest claim we can make 
is ^(0) = h + 0(e). A similar argument applies to higher values of q. 
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5 Analysis of the Full CSP Method 



We now return to the full CSP method and prove Theorem 13.11 Since the 
full CSP method and the one-step CSP method start from the same basis, the 
conditions ()3.14j) and (|4.1()|) are the same for q — 0, 

B\ 0) g = 0. (5.1) 
Therefore, we can choose ip(o) = ^(o) = h Q . 



5.1 Proof of Theorem 13.11 for q = 1 

In this section, we carry out the first iteration of the full CSP method and 
determine the resulting approximation Ke of the slow manifold. We then 
compare Ke and K!P . 

The update quantities C/(o) and L(o) follow from the definition ()3.11j) . 

U m = (AfoV^AfoV ho) = A(o)(A(o)) -1 - (5-2) 

(We recall that we use the same notation U(p) and L(p) for the full matrix and 
the nonzero block.) In particular, Eqs. (|5.2jl and ()4.27|) imply that Z7(o) = £/(o)- 
Next, we update the matrix -B( ). Following Eq. (jSHSJ)? we find 

B (1) = (I-L (0) )(I + U (0) )B (0) . (5.3) 

The upper and lower row blocks of B^ are 

B (i) = B (o) + U (0) B (0)i ( 5 - 4 ) 
Bfa = (J — L(o)C/(o))-B( ) — L(p)B( Q y (5.5) 

Since Pr \ — U(o) = Um\ and -^(o) = VYo)j Eqs. (I4.14|) and (|5.4j) imply that 

so after the first iteration the CSP condition is the same as for the one-step 
method. Therefore, if)m = if)m and, by Theorem 14. X| 

^(i) = h + eh + 0(e 2 ). (5.7) 

This proves Theorem 13.11 for q = 1 . 
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5.2 The Induction Hypothesis 



So far, we have established the identities BL-, = Bhs and BjL = BL from 

which we could conclude that K,f' = icf^ and = K £ • In general, though, 
it is not true that Bis = Bis for higher values of q, as we now demonstrate. 

In the one-step CSP method, Eq. (j4.14|) yields 

B (2) = B (0) + (^(0) + U(l))Bf Q y 

By contrast, in the full CSP method, we obtain from Eq. (j3.13|) 

B (2) = (4 - ^(i)^(o)) -B(o) + (U( ) + U(i) - C/(i)L(o)[/( )) Bf 0) . (5.8) 

The rows of Bhs and B? s are linearly independent, as can be seen from 
Eq. (|3.9jh so the presence of the premultiplier of Bhs in the expression ()5.8|) 

implies that SL 7^ SL. A similar argument shows that Bis 7^ for 
q = 2,3,.... Consequently, the proof of Theorem 13.11 for q = 1 given in 
Section EH] does not generalize to higher values of q. 

The matrix Bis has an important property. Using Eq. (|4. 14)1 . we write 

B \l) = \ P(q-l) B (0) + S (0)> B (0)) = B {0) (^21 B (q-1) B (0) + B {0) > J «) ■ 

Given the induction hypothesis (|4.24j) . we rewrite this expression once more, 



B\ q) = B™ (- g e j Dhj + 0( £ <O, J n j . (5.9) 

Take any y £ K, and let the points Q G Q £ /Ce 9 ™ 1 ^, and Q' G .M e be 

defined by 

Q = (y,i>( q -i)(y,£)), Q = (y,4>( q -i)(y,e)), Q' = (y, h £ (y)). 

The n row vectors of the matrix (—Dh e (y), I n ) form an exact basis for J\fq/A4 £ , 
the space normal to M. £ at Q' . Therefore, by Eq. (|5.9j) . B l ^{Q) is a linear 
combination of the basis vectors of Nqi Ai £ , up to and including terms of 
0(e q ~ 1 ), via the invertible matrix BI&. Hence, the columns of BIJQ) form 
a basis for Mq'M. £ up to and including terms of 0{e q ~ v ). This property of 
Bi s (Q) was central to the proof of Theorem 14.11 We seek to prove a similar 
result for the rows of BlJQ). 

The rows of B^(Q) can be written as linear combinations of the rows of 

B(g)(Q), 

B {q) (Q)=T {q) (y,e)B {q) (Q), (5.10) 
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because B( q )(Q) is invertible (see Eq. (|4.14j) ). In terms of the constituent 
blocks, 

B \l) = T H)B\ q ) + T {q)B(q)' C 5 " 11 ) 

B (g) = T (*)B\>i) + T {%)%)■ ( 5 - 12 ) 

Equation (|5.11|l shows that the requirement that the rows of BhJQ) span 
Mq'M. £ up to and including terms of 0(e q ~ l ) is equivalent to the conditions 

T$(y,e) = 0(1) and invertible, T$(y,e) = 0(e q ). (5.13) 

Assume for the moment that these conditions are satisfied. Then the CSP 
condition (j3.14|) after the gth iteration can be recast as 

r &)(2/> e )£(,)(2/>V'(<i-^ g(y,z,e) = 

or, since T} q \(y,e) is invertible, 

The second term is at least of 0(e q ), by the second assumption in Eq. (j5.13|) . 
so the terms of 0{e^) in Eqs. (|4.10j) and (|5.14|) are equal for j — 0, 1, . . . , q — 1. 

At 0(e q ), the two equations differ by the term (T^^T^Bf^g^y, V>( ? , ), e)- 

Since the 0(1) terms of the two equations agree, it follows that tp( q ,o) = i>(q,o) = 
ho and, therefore, g(y, ip(q,o), e) = 0. Hence, Eqs. (|4.1fi|) and (|5.14jl agree up to 
and including terms of 0(e q ), so Eq. (j5.14j) produces the asymptotic expansion 
of the slow manifold up to and including terms of 0(e q ), by Theorem 14.11 

To complete the proof of Theorem 13-H we need to verify the condi- 
tions ()5.13|) for q = 2,3,..., which we do by induction on q. The induction 
hypothesis is 

T (g) (-,^-i),s) = { £T ^ + 0{£2) }7 +0(£2) J, (5-15) 
^w(-.e) = + (5.16) 

3=0 



5.3 Proof of Theorem [O for q = 2,3,... 

In this section, we carry out the induction step of the proof. We assume that 
Eqs. (|5.15|) and ()5.16|) hold for 0, 1, . . . ,q and prove that they also hold for 
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q + 1. It suffices to establish Eq. f)5.15j) : Eq. ()5.16|) follows immediately from 
Eq. ()5.15|) and our discussion of the CSP condition (j5.14|) . 

Before carrying out the induction step, we derive an update formula for 
T(q). Using Eq. (|5.1()j) with q replaced by q + 1, we obtain 

T (q+1) = B (q+1) A^ q+1 \ (5.17) 

(Here, we used the identity (B^+i))- 1 = A^+V.) Next, we use the update for- 
mulas ()3.13j) and (|4.1j) for and A^ q+1 \ respectively, to rewrite Eq. ([5.17]) . 

T {q+1) = (I- L {q) ) (I + U {q) ) T {q) (I - U (q) ). (5.18) 

Equation (l5~TUJ) also relates A^ q+1 ^ to A^ q+l \ 

= (T^))- 1 . (5.19) 

Taking C = (T( 9 )) 1 in Eq. (j3.19|) . we express in terms of A( 9 ), 



A (9) = T (g) A (g) ( T (9)) 



■(«)■ 



or, equivalently, 



A( 9 ) = r (fl) A (9) (T^)- 1 + (T( g )) 1 (5.20) 



Next, we estimate the blocks of the matrices in Eq. (|5.2(Jj) . The estimate of 
2(g) is given in the induction hypothesis ([5. 15[) : its inverse satisfies a similar 
estimate, 

Also, the induction hypothesis (|5.1fi|) and Theorem 14. II guarantee that ij)^ = 
tp( q ) + 0(e q+1 ), so the expansions of A( g )(y, ip( q ), e) and A^(y,ip^,e) are equal 
up to and including terms of 0(e q ). It follows from Eq. ()4.40|) that 

Taking V = TL) in Lemma IA.21 we conclude from Eq. f[5.15[) that 

_ / 0(. 3 ) O(^) \ 

DT <^-^o( e 2 ) 0(e 3 ) J- (5 ' 23) 
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The desired estimate of A( g ) now follows immediately from Eqs. (|5.15|) . ()5.21|) . 
and dSIH , 

Afi = Ag o) +0( e ), (5-24) 
A(J = ^fc-A^G+O^ 1 ), (5.25) 
Ag) = elAf^+T^All^ + Oie 2 ), (5.26) 
A 22 ) = ^ 2 1) + (9(, 2 ). (5.27) 



The definition (f3~TT|) and Eqs. (j5~23|l and dQHD imply that C/ (?) = 0(e 9 ), with 
the leading-order coefficient given by 

U {q ,g) = (Ago)) 1 Af 2 j(?) = U {m) - TP^. (5.28) 
Furthermore, the definition ()3.11|) and Eqs. ()5.24|) and (j5.26|) imply that 

^)=A 2 i( A (i) _1 = ^)- (5-29) 

Finally, we observe that, to leading order, the blocks of T^(- ,ip^,e) are all 
equal to the corresponding blocks of T^(- ,ip( q -i),e). The latter are given by 
the induction hypothesis (j5.15J) . 

We are now ready to estimate the size of the blocks of T( g+ i)(- , ip^,e). 

The update formula (|5.18|) gives TP +1 * = T}\ + U( q yT?\. According to 
the induction hypothesis, TP = I n + 0(e 2 ) and T? q \ = 0(e). Furthermore, 
U (q) = 0(e q ), by Eq. (pT2%j) . Thus, TP +1) = I n + 0{e 2 ), as desired. 

The update formula dS3H|) also gives TP +1) = TP, - Tf q )U {q) + U {q) Tf 2 - 
U( g )T?\U( g ). According to the induction hypothesis, TP = 0(e q ), TP = 
l n + 0(e 2 ), Tg = 0(e), and T 22 = l m + 0(e 2 ). Furthermore, U (q) = 0{e q ), by 

Eq. (I5~2~%J) . and U {q) = 0(e q ), by Eq. (jOHJ) . Thus, the terms in the formula 
for TP +1 -, are all at least 0(e q ). The same is then true for TP +1 y We will now 
show that TP +1 * is, in fact, at least 0(e q+1 ) by showing that TP +1 q , = 0. To 
leading order, the update formula for TP +1 , is 

T {q+l,q) = T (q,q) ~ U (q,q) + U (l,q)- (5-30) 

Equation ()5.28|) implies that the right member of ()5.30|) vanishes. Therefore, 
T( q+ i t g) = 0, as desired. We emphasize again that the choice of U( q ) is central 
to the working of the CSP method. 



Next, the update formula gives T£ +1) = Tf\-L {q) U iq) T 2 q )-L (q) TP y 

'§ = 0(e) and 



According to the induction hypothesis, T?\ = 0(e) and TP — I n + 0(e 2 
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Furthermore, U {q) = 0{e q ) and L [q) = O(e), by Eqs. lET2S|) and (15311) . Thus, 
the terms in the update formula for T^ +l ^ are all at least 0{e). Hence, T?\ is 
also at least 0(e), as desired. 

Lastly, the update formula (pTTSj) gives T ( 22 +1) = T ( 2 J - L^Tg-Tg^,) + 

L {i) T {qpi.i) ~ L (i) U (g) T (q) + L (i) U (.q) T {qft{q)- According to the induction hy- 
pothesis, T?^ = 7 m + 0(£ 2 ). The remaining terms have already been shown to 
be at least 0(e 2 ). Hence, T 22 +1) = J m + C(e 2 ). 

The proof of Theorem 13.11 is complete. 



6 The Michaelis— Menten— Henri Reaction 

In this section, we apply the CSP method to the Michaelis-Menten-Henri 
(MMH) mechanism of enzyme kinetics to illustrate Theorems 13 . 1 1 and 14 . 1 1 We 
consider the planar system of ODEs for a slow variable s and a fast variable c, 

si = e(-s + (s + K- X)c), (6.1) 
c' = s-(s + k)c. (6.2) 

The parameters satisfy the inequalities < e 1 and k > A > 0. Only 
nonnegative values of s and c are relevant. The system of Eqs. ()6.1)) - (j6.2|) is of 
the form ()2.1j) - (j2.2|) with m — 1, n — 1, y — s, z — c, g\ — — s + (s + k — A)c, 
and g 2 = s — (s + k)c. 

In the limit as e j 0, the dynamics of the MMH equations are confined to 
the reduced slow manifold 

M = {(c,s):c=^—,s>0}. (6.3) 
s + n 

The manifold M.q is normally hyperbolic, so according to Theorem 12.11 there 
exists, for all sufficiently small e, a slow manifold M. £ that is 0(e) close to Aio 
on any compact set. Moreover, A4 e is the graph of a function h £ , 

M £ = {(c,s):c = h £ (s),s>0}, (6.4) 

and h e admits an asymptotic expansion h e = h + ehi + e 2 h 2 + • ■ • . The 
coefficients are found from the invariance equation, 

s-(s + K)h £ (s) = £h' £ (s)(-s + (s + k- X)h £ (s)). (6.5) 

The first few coefficients are 

... s , . . k\s , , . kXs(2kX — 3Xs — ks — k 2 ) 

ho{s) = —.> l{s) = 2{) = M • 

(6.6) 
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6.1 Application of the One-Step CSP Method 



Both the one-step and two-step CSP methods start from the same initial basis. 
We choose the stoichiometric vectors as the basis vectors, so 

•4<o, = (4°>, 4°») = ( ? i y % = c so = c ? i y (6.7) 



The CSP condition Bh^g = is satisfied if c = h (s), so the CSP manifolds 
£e and /Ci '* coincide with Ado- With this choice of initial basis, we have 



A {0) = B {0) (Dg)A^ = ( 



-(s + k) -(c-1) 
e(s + « — A) e(c — 1) 



(6.1 



First iteration. At any point (s, c) , we have 

^ (1) = (? _M. A.» = (t J )■ < 8 - 9 > 



On IC^\ these expressions reduce to 

1 



= ( i ^ ) • = ( "r n ) ■ ( 61 °) 



The CSP condition, 



~1 , s k( — S + (s + K — X)c) 

B l {l) g = 3-(s + k)c- ^— -- } -L = 0, (6.11) 



is satisfied if 



s k\s 9 k 2 Xs(s + k — X) „, ,. 



s + k (s + k) 4 (s + k) 



Comparing this result with Eq. (|6.6|) . we see that the asymptotic expansions 
of Ks and M. E coincide up to and including 0(e) terms, in accordance with 
Theorem 14. II for q = 1; however, the 0(£ 2 ) terms differ at this stage. 



Second iteration. The blocks of Am are 

Aft = - (s + K) +£ (g + K : A)(c - 1) , (6.13) 



S + K 

xia - „ | Jc-1)[A(c-1)-(-s + (s + k-A)c)] 

Aft = e(s + « — A), Ag=3^. (6.10) 

S + K 
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On K!i\ the blocks reduce to 

111 / \ k(s + K — \) 2 K\s(s + K — A) 

a ;;» = -(.+«>-« (S+K)2 +e (.+«)» ■ (616) 

ri2 _ kA(k-2s) 2 kAs(2k(s + k - 2A) + As) 

A 21 = £(s + K _ A); A 22 £ ^^ +£ 2^^_ (61g) 

The second update is 

A? = (J), (6.19) 

= ( g J + £ I «A(<f-3s) 

+ £ 2 ( kA[k(5s-k)(s+k-A)+As(s-2k)] ) +0(s 3 ), (6.20) 

2 / kA[k(5s — «)(s + k — A) + As(s — 2k)] 

+ £ V (^F 

+ 0(e 3 ), (6.21) 
Bf 2) = (1,0). (6.22) 

The CSP condition 

51 / . \ k(-s + (s + k- X)c) 
B\ 2) g = s-(s + k)c-s 

(s + /t) 5 v ' 



is satisfied if 



(6.23) 



s k\s kXs(2k\ — 3Xs — ks — k 2 ) „. 

c = —. + "v^ +e (.+«)' * +<*•)•(•■"> 

Comparing this result with Eq. (|6.6|) . we see that the asymptotic expansions 
of Ke and *M e coincide up to and including 0(e 2 ) terms, in accordance with 
Theorem 14. II for q = 2. 
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6.2 Application of the Full CSP Method 



First iteration. At any point (s,c), we have 



(i) 



On K,f\ these quantities reduce to 







S + K — X 



1 / S + K 

S + K I 



C-l } 7 

S + K 



A? 



c-l J , 

S + K 



B 2 {1] = (1, 0) + e 



s + k — A / c — 1 



S + K V S + K 



(6.25) 
(6.26) 



B \i) 



(i) 







S + K — A ( 1 



K I 5 



(1) 



(S + K) 2: 



S + K \ ( S + K ) 



1), Bfa = (l,0) + e 



A 



S + K — X 

S + K 



K I ) 

(s+/t) 2 



K 



(6.27) 



:S + k) 2 ' 



The matrix relating to its one-step counterpart is 



T, 



(i) 



1 

S + K — A 



S+K 

(0) 



(6.28) 



(6.29) 



so T(i) is indeed of the form (|5.15|) on K\ . 

Equations ()6.10|) and ()6.28|) imply that Bm = Bm , so the CSP condition 
yields ipm = ipm. Thus, after one iteration, the full CSP method also finds 
the expansion of A4 £ up to and including 0(e) terms. 



Second iteration. The blocks of Am are 



A 11 

A (i) 



A 12 

A (l) 



A 21 
A (l) 



Is + K) + s 



[s + K — A) 



+ e 
s 

S + K 

e 2 



S + K 
2 V C — l)(s + K — A) 



(c - 1) + (c - -f-) 



— C + £ 



(s + k) 3 
c-l 



S + K 

-A(c - 1) + (-s + (s + k - A)c)] ,(6.30) 



[A(c - 1) - (s + (s + k - A)c)] , (6.31) 



fs + k) 2 



v s + k) 2 
c - l)(s + « - A)(s + k - 2A) 



+A(-s + (s + k - A)c) + (s + k - A) 2 c 



S + K 



(6.32) 
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A22 

A (l) - 



A(c-1) + (s + k-A)( 



S + K 

_ 2 (c-1)(s + k- A) 



c) 



S + K 

[A(c-l)-(-s + (s + «-A)c)] 



with remainders of 0(£ 3 ). On /C^, the blocks reduce to 



A 11 

A (i) 



A 12 

A (l) 



A 21 

A (l) 



a 22 
A (l) 



= e 



, s k(s + k — X) 9 kX(s + k — A)(3s — k) 

-(s + k) - e + — -. ^ 

v ' (s + k) 2 (s + k) 5 

kX(k — 2s) 2 kXs(2k(s + k — 2A) + As) 

(s + K ) 4 +5 (s~+~kJ ' 

2 k(s + k - A)(s + re - 2A) + A 2 s 



— £ 



— £■ 



kX 



s + re) 2 



— e 



(s + k) 3 

2 kX((2s — re)(s + re — A) — As) 
(s + re) 5 



with errors of 0(e 3 ). The result of the second iteration is 



A (2) 
^12 



s + re — A 2 re(s + re — 2A)(s + re — A) + A 2 s 



A (2) _ 
^21 _ 



1 + e 
1 - e 



S + K (s + re) 4 

2 reA(2s — re)(s + re — A) 



(s + re) 6 
re(s + K — A) 
(s + re) 3 

2 (s + re - A)(re 2 (s + re - 2A) + reAs) + kX 2 s 



(s + k) 



A 



(2) 
22 



K 



(s + k) 2 

+ € 



+ 6 



reA(re — 3s) 
(s + re) 5 



2 k 2 X(7s - 2k)(s + re - A) + reA 2 s(s - 2k) 



(s + k) 



R 11 



(s + k) 



k kX(k — 3s) 

— e- 



(s + k) 1 



B 



12 

(2) 



R 21 



l + e 
l-e 

+ € 



,re 2 A(7s - 2k)(s + k - A) + reA 2 s(s - 2k) 
' (s + re) 8 ' 

2 reA(2s — re)(s + re — A) 
(s + re) 6 ' 
re(s + re — A) 



(s + re) 3 

2K s + k — X)(k 2 (s + k — 2X) + reAs) + reA 2 s 



[s + kY 
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22 s + k-X 2 k(s + k- A) (a + k-2A) + A 2 s 

5 ( 2 ) = £ ^T^" £ fT^ ' (6 - 45) 

up to and including terms of 0(e 2 ). Also, on K?i\ 
u ^ 2 kA(s + k- A)(2s-k) 

21 (s + k — A) 2 K ( S + K — ^)( S + K — 2A) + A 2 s 



T'22 -i 2 

(2) — 1 ~ £ 



k\(2s — k)(s + k — A) 



,5 + K) a 

with remainders of 0(e 3 ). Thus, T( 2 ) is indeed of the form (|5.15|) on jdP ■ 
The CSP condition 

k{— s + (s + k — A)c) 



5^.9 = s - (s + k)c - e 
+ e 2 n\ 



(3s — k)(— s + (s + k — A)c) 



is satisfied if 



(s + k) 5 

(2s — k)(s + k — A)(s — (s + k)c) 

+ (7+^p 

+ 0(e 3 ) 

0, (6.46) 



s kAs 9 kAs(2kA — 3As — ks — k 2 ) „. , 
+ ^7 u + e " 7 w + 0(e 3 ). (6-47) 



s + k (s + k) 4 (s + k) 7 

Therefore, after two iterations, the full CSP method finds the expansion of 
M. £ up to and including 0(e 2 ) terms. 



6.3 The Second Step and the Fast Fibers of M £ 

The preceding analysis of the full CSP method shows that, at the qth iteration, 
the second step alters only the terms of 0(e q+l ), leaving the terms of 0(1) 
through 0(e q ) invariant. Here, we observe that the second step also plays a 
constructive role for the dynamics near the slow manifold. As can be seen 
in the MMH example, the second step yields the asymptotic expansions of 
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the tangent spaces of the fast fibers at their basepoints up to and including 
terms of 0(e q+1 ), at least for q = 0,1, and 2. This additional information is 
contained in the columns of We remark here that this property is not 

shared by the one-step CSP method, since the columns of remain tangent 
to the fast fibers at their basepoints only to leading order after each iteration. 
Details about the fast fibers and their tangent spaces will be presented in a 
future publication. 



7 Relation between CSPM and ILDM 

The CSP iteration procedure is designed to diagonalize the Lie bracket [-,<?]. 
At each iteration, the then-current basis is updated in such a way that [ • , g] is 
block-diagonalized to the next-higher order in e. Thus, each iteration improves 
the quality of the basis of the orthogonal complement of the tangent space. 
The CSPM is defined as the locus of points where the vector field is orthogonal 
to that orthogonal complement. 

The ILDM method works, instead, with the Jacobian, Dg, of Eqs. (j2.3)l - 
(12.4)1 . A Schur decomposition transforms Dg into upper triangular form, 

Dg = QNQ>, N=( N Q S N ^ ) , (7.1) 

where Q = (Q s Qf) is unitary. The eigenvalues of Dg appear on the diagonal 
of iV in descending order of their real parts, from least negative in the upper 
left to most negative in the lower right. The first m Schur vectors (the columns 
of Qs) form an orthogonal basis of the slow subspace and the remaining n Schur 
vectors (the columns of Qf) an orthogonal basis of the orthogonal complement 
of the slow subspace. The vector field g is entirely in the slow subspace if it is 
orthogonal to this orthogonal complement — that is, if 

Q' f 9 = 0. (7.2) 

This equation defines the ILDM; see [TTJ Section 3]. 

As we showed in jllj . the ILDM is only a first-order approximation to M. £ . 
The error is always 0(e 2 ) unless Aio is linear. The error can be traced back 
to the choice of the operator. The tangent space is a left-invariant subspace 
of the Jacobian only to leading order, so putting Dg in upper triangular form 
yields the orthogonal complement only to leading order. Since the linearized 
system is only an approximation of the original ODEs ()2.1j) - ([2.2p . this choice 
does not produce an exact result unless g is linear. The success of the CSP 
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method in approximating the slow manifold is due to the fact that the ODEs 
for the amplitudes / are equivalent to the ODEs (|2.1|) - (j2.2j) . That is, the full 
nonlinearity is retained. 

The time-derivative term in the definition ()3.tij) must be included in the 
evaluation of A; otherwise, the accuracy of the CSP method is compromised. 
In fact, such an omission results in implementing the ILDM rather than the 
CSP method, which may be seen as follows. With our initial choice of a point- 
independent basis A(°\ the matrix A( ) is similar to Dg; see Eq. (j3.10j) . The 
omission of the term (dB( q )/dt)A^ in the calculation of A( 9 ), for q = 1,2,..., 
would lead to the formula A( g ) = (I + P^)B^(Dg)A^°'(I — P( q )), which would 
imply that A( g ) is similar to Dg. Therefore, the one-step CSP method would 
put Dg, rather than A, in lower-triangular form, just like the ILDM method. 
After the second iteration, one would make an error (proportional to the cur- 
vature of M.q) at 0(e 2 ), which subsequent iterations would not remove. The 
MMH example in Section |U] illustrates these observations. 
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A Auxiliary Lemmas 



Lemma A.l The quantity Dh q is given by the formula 



Dh n 



9-1 



9-1 



(D y g 2 ) q + J2 P*ftO,-< Dht - J2 D " h 



£9l,q-£-l 



9-1 



i=0 

q—1 q—l—i 



e=o 



j2(Dh e )(D ygi ) q ^ e -j2 E Dh * D 

z9i) q -i-i-eDhi 



i=0 1=0 



(A.3) 
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Proof. The coefficient h q is found from the 0{e q ) terms in the invariance 
equation (|2.10p . 

q-l 

g 2>q = J2( Dh i)9i >g -e-i- (A-4) 

1=0 

Taking the total derivative with respect to y of both sides of (jA.4|) . we find 



d q ~ l 9-1 d 

-f92, q = ^2(D 2 h e )g hq _ e _ 1 + ^2(Dh l ,)—g hq _ l ,_ 1 . (A.5) 
^ 1=0 e=o y 

The operations of taking the total derivative with respect to y and expand- 
ing with respect to e commute, because the Fenichel theory guarantees C r 
smoothness in e and y for each r. Therefore, 

d ( dg, x 



92, q = (^) ={D y g 2 ) q + y £(D,g 2 )^(Dh i ), (A.6) 
^ y ' q i=o 

9i,q-i-i = (^) =( D v9i) q -i-i+ E (^i) ? -i-^(^)-(A.7) 



d ( dg 

d y \»*/ q -i-i - i=0 

Substituting Eqs. ()A.6|) and (jA.7j) into Eq. (|A.5|) . we obtain 

q q-i q-i 

(D y g 2 ) q + J2^ D ^) q - l Dh l = g^^^H + g^KD^),^ 

i=0 i=0 t=0 

q-l q-l-l 

+ E E {Dh^D.g^^Dhi). (A.8) 

t=0 i=0 

Separating the i — q term in the sum of the left member, changing the order 
of summation in the last sum of the right member, and solving for Dh q , we 
obtain Eq. (|A~3|) . I 



Lemma A. 2 Let V be a matrix-valued function ofy, z, and e that, together 
with its first-order derivatives, is smooth and 0(1) as e \ 0. If z — ip^(y,e) 
and 

q q 
V(-,i> ig) ,E)=J2 £iV t + °( £q+1 )> 9i(-,i>( q) ,e) = J2 £e 9i,t+0(£ 9+1 ), (A.9) 

i=0 1=0 

then, 

% (• . = E ^ E ^9i,i-i + 0(e^). (A.10) 

i=0 t=o " 
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Proof. A direct computation gives 

^ = (DV)g = e(D y V) 9l + (D z V)g 2 , (A.ll) 

where all the terms are evaluated at (y,ipf q \(y,e),e). Since ipt q ^ approximates 
the slow manifold up to and including 0(e q ) terms, 

9i(-,i>( q) ,e) = 9i(-,K,e) + 0(e q+1 ), (A.12) 
g 2 {- , V (g) , e) = g 2 (- , h £ , e) + (A.13) 

and also 

DiP (q) = Dh £ + 0(e q+1 ). (A.14) 
Using Eqs. (ETTUJ) . (fA~T2ll. and (TO!) , we rewrite Eq. (lA~T3l) as 

#,(• ,V (g) ,e) = e{D^ {q) ) gi {- ,V (g) ,e) + (A.15) 

Equation (jA.15|) is an equation for /Q 9 . We recast it so the right member 
involves a total derivative with respect to y, 

(DV)g = e {D y V + D z VD^ [q) ) 9l + O^ 1 ) = e^- 9l + 0{e q+l ) (A.16) 

or, expanding in powers of e, 

{DV)g = E^ +1 E(^) 9^ + 0(s^). (A.17) 

The operations of taking the total derivative with respect to y and expanding 
with respect to e commute, so (dV/dy)e = dVe/dy and Eq. (|A.10|) follows. I 
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